Progression of Phase Behavior for a Sequence of Model Core-Softened Potentials 



in 

O 

o 

(N 

» 

o 
O 

(N 



O 



D QuiglejQ and MIJ Probert 

University of York 
Heslington 
York 
United Kingdom 
YO10 5DD 
(Dated: February 2, 2008) 

A series of phase diagrams is obtained in three dimensions for a smooth pair-potential with an 
outer well and a repulsive inner shoulder. Condensed phase boundaries are located using free energy 
calculations. Liquid-vapour equilibria are obtained with multi-canonical methods. As the depth of 
the outer well is increased, a simple-hexagonal to close packed transition appears in the solid leading 
to a discontinuity in the slope of the melting curve. For deeper wells the simple hexagonal melting 
temperature exhibits a maximum with respect to pressure. The onset of the predicted metastable 
isostructural transition is also studied. 

PACS numbers: 61.20.Ja, 61.66.Bi, 64.70.Dv, 64.70.Fx 
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I. INTRODUCTION 

Since the work of Stell and Hemmer [lj it has been 
appreciated that a simple model pair-potential with two 
repulsive regions of differing strength is capable of gen- 
crating a third fluid phase. Recent experimental and 
theoretical evidence for liquid-liquid phase transitions 
(LLPT) in elemental melts 0, H H has revived in- 
terest in these models as a mechanism for reproduction 
and study of the gener al phenomenon. In pa rticular, it 
has been suggested bv iMishima and Stanley! [(J that so 
called 'core-softened' potentials are capable of generating 
a first-order phase transition within a supercooled liquid. 
LLPTs of this kind are observed in many models of water 
011 H [Hill 113 and may be related to the celebrated 
density maximum at 4°C. The LLPT is considered an 
extension of the transition between two amorphous ice 
phases into the supercooled liquid regime, ending in a 
critical point below the thermodynamic melting temper- 
ature. 

A continuous interpretation of a Stell-Hcmmcr like po- 
tential can be reproduced by the function 
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which for appropriate choices of A, r n and w closely 
resembles t he form suggested by IMishima and Stanley! 
IScala et all [l3j have presented a clear argument that 
potentials of this form posses an isostructural transition 
between either two solid or two liquid phases of different 
density. Previous simulation studies have employed both 
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FIG. 1: Core softened pair-potentials. The discrete form is 
shown as a dashed line, along with a corresponding smooth 
form (solid line) constructed as in equation The Lennard- 
Jones potential (dotted line) is shown for comparison. 



the smooth form given by equation^and a discrete piece- 
wise interpretation. Both are shown in Fig. ^ I n two di- 
mensions the observation of a water-like density anomaly 
has been reported in molecular dynamics simulations of 
both the smooth and discrete forms 0, 0] • Other sim- 
ulations have employed advanced Monte-Carlo methods 
[ltij . These indicate that rather than being related to a 
metastable LLPT, the density anomaly is a consequence 
of cluster formation during quasi-continuous freezing to 
a solid of lower density than the surrounding liquid. 

Simulations in three dimensions have been largely 
limited to the discrete potential, and indicate a sec- 
ond critical point can occur for certain parameterisa- 
tions 0, m> m> H0| > without a density anomaly. This 
metastable second critical point can lie at higher or lower 
temperature than the liquid-gas critical point depending 
upon the specific parameterisation employed. 

Previous simulations of the smooth potential in three 
dimension have been restricted to two limited studies 
pll |23 |. In a recent rapid communication [2j|, we 
mapped the low pressure portion of the three-dimensional 
phase diagram for A — e. A simple hexagonal (sh) 
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so lid phase, identifi ed using the meta-dynamics method 
of iMartonak et ah! |24|. was shown to remain thermody- 
namically stable in this regime. In this paper we report 
on calculations for a range of A values, studying the emer- 
gence of the sh phase. The sh to close packed transition 
at high pressures is studied. In addition we follow the 
behaviour of the melting and liquid-vapour critical tem- 
peratures as a function of the outer well depth. The 
possibility of isostructural phase transitions is also inves- 
tigated. 

Measured quantities are presented here in the usual di- 
mensionless reduced units. Energies are quoted as mul- 
tiples of the inner-well depth e, lengths as multiples of 
a. Reduced temperature T* is calculated as fc^T/e, 
with pressure P* = Pa 3 /e. Time is measured in units 
t* = (m/e) 1/,2 (T where m is the atomic mass. For all 
studies reported here the parameters = 1.44<t and 
w = 41.22<7~ 2 are used in equation^ The pair-potential 
is truncated at 2.5a and /orce-shifted such that no dis- 
continuities occur in the dynamics. 

The remainder of this paper is organised as follows. 
In section ITU we explore the zero temperature phase be- 
haviour as a function of the outer well depth A and choose 
values for study at finite temperature. In section UTTl we 
describe the methods employed for our finite temperature 
studies. Results are presented in section llVl Finally, we 
summarise our findings. 



II. ZERO TEMPERATURE BEHAVIOUR 

Much insight into the phase behaviour of a model sub- 
stance can be extracted from zero temperature enthalpy 
and volume characteristics. The zero temperature en- 
ergy is obtained by performing a conjugate-gradient (CG) 
minimisation of the potential energy U with respect to 
both the cell vectors and the fractional atomic coordi- 
nates. Our identification of candidate structures was re- 
ported in Ref. 23]. The energy under zero pressure for 
energetically relevant structures is plotted as a function 
of the outer well depth A in Fig. [21 In order of decreas- 
ing density, these are the high density face centred and 
hexagonal close packed structures (hd-fcc and hd-hcp) 
associated with the unmodified Lennard-Jones potential, 
simple hexagonal (sh), simple cubic (sc) and a second set 
of lower density close packed structures (ld-fcc and ld- 
hcp). The final two structures are stabilised by the outer 
well only for A > 0.7 e. 

The energetic favourability of open structures on in- 
creasing A can be understood in terms of neighbour dis- 
tances. The sh primitive cell is defined by a = b = 1.059<t 
and c = 1.016ct with a = /3 = 90°,7 = 120°. The nearest 
and second nearest neighbours lie at a distance c and a 
respectively, close to the position of the Lennard -Jones 
minimum, while the third nearest lie at yj (a 2 + c 2 ) which 
is very close to Tq. Both energy minima in the pair po- 
tential are therefore utilised. In the sc structure, nearest 
neighbours lie at approximately 1.12a. The favourability 
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FIG. 2: Optimised energy at zero temperature and pressure 
for various structures as a function of the outer well depth 
parameter. The hep and fee structures are near-degenerate 
as indicated in the inset. 



of sc over the high density close packed structures stems 
from second nearest neighbours which lie at times 
this distance, which is close to tq. Each atom has six 
neighbours close to the Lennard-Jones minimum plus six 
close to r . This compares to eight and twelve in the sh 
case which is therefore lower in energy. 

Both pairs of close-packed structures decrease in en- 
ergy on increasing A. The lower density structures con- 
tain more neighbours close to ro and hence the decrease 
is faster in these cases. 



III. METHODS 
A. Liquid- Vapour Equilibria 

Liquid-vapour equilibria are traced in the direction of 
decreasing temperature from near the critical point. An 
initial estimate of this critical point is obtained using 
molecular dynamics simulations in the canonical ensem- 
ble. Fluid isotherms are traced in the pressure-density 
plane at temperatures increasing in steps of AT* = 0.08. 
The temperature of highest isotherm exhibiting the clas- 
sic hysteresis associated with the liquid- vapour transition 
is taken as an estimate of the critical temperature. Typ- 
ically between 4 and 8 isotherms are required for this 
process. Each isotherm is computed from 30 simulations 
over the density range p* = 0.1 to 0.8. Simulation time 
at each density is t* — 400 in both equilibration and 
production periods. 

Tracing of the liquid-vapour coexistence curve pro- 
ceeds using Monte-Carlo simulations in the grand- 
canonical ensemble (GCE). These combine multi canon- 
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ical samplin g 1251 with re-weighting of particle number 
histograms |2g. O ur metho d closely follows that de- 
scribed in detail by IWildinel |27j . A starting point on 
the liquid-vapour coexistence curve is obtained from the 
above estimate of critical temperature by fine tuning of 
the chemical potential fi* until a bimodal distribution is 
measured in the particle number N, with equal area in 
each peak. This is then re-weighted to a lower temper- 
ature, providing a biasing function for multi-canonical 
sampling. The process repeats, proceeding down the co- 
existence curve in steps of AT* = —0.008. At each step a 
multi-canonical GCE simulation of 500 000 Monte-Carlo 
cycles is sampled at equilibrium to produce the parti- 
cle number histogram. A cubic simulation cell of side 
L = 7.13cr is employed in all cases. Typically twenty 
steps are employed, tracing the curve over a tempera- 
ture range of approximately AT* = 0.15. The resulting 
data is used to estimate the parameters of the liquid-gas 
critical point. 

At lower temperatures the liquid-vapour curve is 
traced by numerical integration of the Clausius- 
Clapeyron equation (see below). 



B. Solid-Solid Transitions 

These are located via free energy calculation of both 
phases along an appropriately chosen isotherm. Free en- 
ergies are calculated using the Einstein crystal method of 
Frenkel and Ladd [28| . Here the free energy derivative is 
integrated by sampling typically 20 — 50 points along the 
path connecting the core-softened solid to the reference 
harmonic crystal, using canonical ensemble Langevin dy- 
namics simulations. Each sample employs a total simu- 
lation time t* = 3.6 at equilibrium. The appropriate cell 
in which to perform this thermodynamic integration is 
first identified by employing constant pressure Langevin 
dynamics simulations |29| of typical duration t* = 100 
after equilibration. A fully flexible simulation cell is used 
in all cases. System sizes of N = 256 and 392 use used 
for fee, sh structures respectively. 

Two sources of system size dependence arise in these 
calculations. The first arises in the thermodynamic inte- 
gration itself. This has been largely compensated for with 
repeat integrations at larger N. The free energy is then 
ex trapolated to the thermodynamic limit via the method 
of iPolson et alJ [30| . The second arises from finite-size 
effects in the identified density at a given pressure. A 
selection of repeat free energy calculations at densities 
obtained from larger system sizes have shown this to be 
negligible. In both cases the additional system sizes used 
are N = 500 and 864 for fee, N = 640 and 864 for sh. 

Free energy calculations of hep and sc structures are 
also reported below, but are not used in locating phase 
boundaries. No finite size analysis has been performed in 
these cases. System sizes used are N = 384 for hep, and 
N = 343 for sc. 



C. Melting Curves 

Melting curves are identified by two independent meth- 
ods in this work. The first is based on direct simulation 
of phase coexistence in the NPT and NPH ensembles 
|31j. The NPT ensemble method is used here to study 
melting of fee structures. A cell is prepared with 500 
atoms in each phase and simulated at a specified tem- 
perature and pressure for t* — 300. During this time the 
system transforms to a pure solid or liquid state. A series 
of these simulations allow the melting temperature to be 
accurately bracketed. 

In contrast, the NPH method requires a single sim- 
ulation to locate the melting temperature at a specified 
pressure ■ Provided the initial enthalpy of the system 
lies close to that at the melting temperature, fractions 
of the system will melt or solidify, shifting portions of 
the conserved enthalpy between the two phases. The 
thermodynamic melting temperature is obtained at equi- 
librium. We note that NPH simulations employing an 
Andersen-Hoover [33ll3^ | or (as in this work ) a Martyna- 
Tobias-Klein [3^| barostat, the enthalpy is not exactly 
conserved, but is constant to within fluctuations in the 
fictitious kinetic energy associated with the cell dynam- 
ics. This pseudo-NPH approach is therefore only useful 
in cases where this fluctuation can be kept within the 
latent heat per atom associated with the melting transi- 
tion. This criteria has only achievable for melting of the 
simple hexagonal structure. In this case a simulation cell 
with 332 atoms in each phase is prepared and simulated 
as described in Ref. [23|. 

Melting temperatures for each potential are also com- 
puted from solid and liquid free energies along an appro- 
priately chosen isobar. Solid free energies are computed 
as above with finite-size corrections applied in all cases, 
up to the maximum temperature of mechanical stability. 
Liquid free energies are computed via thermodynamic in- 
tegration from a reference point of known chemical po- 
tential. This is obtained as follows. 

A point is chosen within the liquid region of the T — [i 
plane. The choice is aided by the data obtained during 
the multi-canonical sampling, and by the solid free en- 
ergy calculations which yield the chemical potential of 
the superheated solid. At the chosen reference point a 
further GCE Monte-Carlo simulation is conducted to ob- 
tain density and pressure information. The Helmholtz 
free energy at the reference point can then be calculated 
to within the statistical uncertainty inherent in this sim- 
ulation. These GCE simulations employ 500 000 MC 
cycles in a cubic simulation cell of side 7.13er. 

Liquid free energies along the isobar of interest are ob- 
tained using thermodynamic integration to the density 
(identified using constant pressure Langevin dynamics 
simulations as for the solid) and temperature of inter- 
est. Each integration samples the relevant free energy 
derivative at 10 points along an isotherm and isochore 
using canonical ensemble Langevin dynamics simulations 
of duration t* = 50. A system size of N = 350 is used. 
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Errors in the thermodynamic integration have been es- 
timated by evaluating the change in free energy around 
various closed loops in the T — p plane. The system 
size dependence has been estimated by repeating a se- 
lection of integrations with N = 600, and by employing 
larger simulation cells in GCE simulations. The domi- 
nant source of error is identified as statistical uncertainty 
in the liquid reference point. This has been controlled by 
employing suitably long GCE simulations. 



D. Clausius-Clapeyron Integration 

The free energy calculations described above locate a 
single point on a phase boundary. To trace the remain- 
der of a phase coexistence curve, the Gibbs-Duhem in- 
tegration method j^| is employed. Beginning from the 
identified single point, the Clausius-Clapeyron equation 
is evaluated over two (one for each phase) constant pres- 
sure Langevin dynamics simulations. This is then in- 
tegrated with the fourth-order Runge-Kutta method to 
trace the phase boundary in the P—T plane. Simulations 
at each point are typically of duration t* = 100 to 150. 
System sizes are N = 500 for fee phases and N = 640 for 
sh phases. In the case of melting curves the number of 
atoms in the liquid matches that in the solid. 

For tracing liquid-vapour coexistence curves a system 
size of N = 500 is used for both phases. An initial point 
for the series is taken from a fit to temperature-pressure 
data unfolded from the multi-canonical GCE simulations. 

Integration error is estimated by reversing each scries, 
integrating back toward the initial point. In all cases this 
error is much smaller than the uncertainty of the initial 
point. 




(a) 
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(b) 



IV. RESULTS 

A. Phase Behaviour for A = e/4 

At A = e/4 the outer Gaussian well represents a 
small perturbation on the Lennard- Jones potential. No 
significant change in phase behaviour is therefore ex- 
pected. Following the procedure in section fill Al an ini- 
tial bimodal particle number histogram was obtained at 
T* = 1.100 with a chemical potential of p* = -9.267. 
The liquid-vapour curve was traced to a temperature of 
T* = 0. 925 i n 21 steps. Example histograms are shown 
in Fig 



FIG. 3: Liquid-gas transition in the A — e/4 potential. Sam- 
ple histograms obtained during multi-canonical GCE simula- 
tions are shown in (a). At lower temperatures the transition 
is traced using Gibbs-Duhem integration. The join between 
the two methods in the P — T plane is shown in (b). 



This process identifies the critical point at T* — 1.108 ± 
0.003, p* = 0.302 ± 0.003. Note that the finite-size de- 
pendence of these parameters has not been investigated. 

Below T* = 0.925 the continuation of the liquid- vapour 
coexistence curve has been traced with Gibbs-Duhem in- 



3§]] These yield information on the density of tegration. The continuation is shown in Fig. [ggg 



the two phases, allowing a fit to the scaling law 

p\ -p* g <x (T* - T* c f 
and to the law of rectilinear diameters, 



Pi 



p* + A(T* - T*). 



(2) 



(3) 



Free energy calculations have been performed at tem- 
peratures in the range T* = 0.083 to 0.500 in steps of 
0.042 along the zero pressure isobar. At no temperature 
in this range were the sc and sh structures mechanically 
stable. As the fee structure is of the highest density, we 
can conclude that as with the Lennard- Jones system, the 
solid remains close packed over the entire positive pres- 
sure range. 
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FIG. 4: Melting curve for the potential with A = e/4 calcu- 
lated using the two-phase NPT coexistence method described 
in section IHI CI The dashed line is the melting curve com- 
puted from free energy calculations and traced with Gibbs- 
Duhem integration. 



Results of two-phase NPT melting calculations close 
to the melting line are shown in Fig. 0] 80 simulations 
were conducted in total. Each sampled point is visually 
identified as either solid or liquid. Simulations which 
could not be identified as a pure phase within the simu- 
lation time (i.e. those close to the melting line) are not 
plotted. 

The melting temperature at a pressure of P* — 0.047 
has been calculated accurately by free energy calculation. 
For the solid phase, calculations were performed a tem- 
peratures up to T* — 0.500 in steps of 0.042. A reference 
point for liquid free energy calculations was taken at T* = 
1.000 and total chemical potential p,* = -10.526. GCE 
Monte-Carlo simulation allows the Helmholtz free energy 
at this point to be identified as /* = —11.080 ±0.003 per 
atom. Six free energies along the P* = 0.047 isobar have 
been computed by thermodynamic integration from this 
reference point, at temperatures of T* = 0.375 to 0.583. 
Interpolation to the intersection reveals a melting tem- 
perature of T* — 0.471. The total error on this melting 
temperature is estimated at less than ±0.008, dominated 
by statistical uncertainty in the pressure and density of 
the liquid reference point. 

This melting point has been used to begin Gibbs- 
Duhem integration in the direction of increasing temper- 
ature. The resulting series is shown in Fig. 0|and in Fig. 
along with the liquid- vapour coexistence line. 



B. Phase Behaviour for A — e/2 

For this value of A we are interested in determining 
if the sh structure remains thermodynamically stable at 
appreciable temperature, and the location of the possible 
sh-fee transition. 

At a temperature of T* = 1.279 a bimodal his- 
togram on the liquid-vapour coexistence curve was ob- 
tained at fi* = —10.782. Subsequent histograms ob- 




1.2 



0.8 - 



0.6 



- 1 : 


i 


_J \ liquid • • 






fee 


liquid • ; , : 


solid 


— + vapour '. : ; 




1 ) : l\ 


i 



0.5 



1.5 



FIG. 5: Phase diagram of the A = e/4 potential (inset). 
The temperature-pressure and density-temperature projec- 
tions are shown. Both forward and reverse Gibbs-Duhem se- 
ries are plotted. The two are indistinguishable on this scale. 
The low pressure region and triple point have not been studied 
in detail. Other than shifts in the melting and critical temper- 
atures, no interesting phase behaviour over the Lennard- Jones 
case is observed. 



taincd along the transition identify the critical point at 
T* = 1.293 ± 0.006, p* = 0.284 ± 0.006. Below temper- 
atures of T* — 1.142 the liquid- vapour curve was traced 
with Gibbs-Duhem integration. 

The resulting free energies are shown in Fig. [BJ These 
are uncorrected for finite-size effects. The sc structure is 
now mechanically stable at finite temperature, but not 
beyond temperatures of T* = 0.25. Although at zero 
temperature the energy of the sh structure is slightly 
less than that of the fee, at no finite temperature re- 
alisable in a simulation is this the case. As temperature 
increases the separation between the fee and sh structures 
increases. There must therefore be a transition between 
these two structures. The free energy difference is how- 
ever too low to resolve with the methods employed here. 

A total of 63 two-phase NPT simulations locate the 
melting temperature between T* = 0.3 and T* = 0.4 
at pressures P* < 0.7. An accurate point on the melt- 
ing curve was sought along the P* = 0.047 isobar. For 
the solid phase, free energy calculations were performed 
for temperatures in the range 0.208 to 0.5 in steps of 
0.042. For liquid free energy calculations, a reference 
point at T* = 1.167 with /i* = -9.474 was taken with a 
Helmholtz free energy per atom of /* = —12.75 ± 0.01. 
Thermodynamic integration to temperatures between 
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FIG. 6: Gibbs free energy per atom as a function of tempera- 
ture along the zero pressure isobar for the A — e/2 potential. 
Three structures are shown. The sc structure is mechanically 
unstable beyond T* = 0.25. The fee and sh structures are 
indistinguishable over much of the temperature range. 
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FIG. 8: Phase diagram for the potential with A = 0.55e in 
the region of the sh-fcc-liquid triple point. The three Gibbs- 
Duhem series form a triangle at their intersection indicating 
the area in which the triple point is located. 
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FIG. 7: Phase diagram of the A = e/2 potential (in- 
set). The temperature- pressure projection and density- 
temperature projection are shown. This phase diagram repre- 
sents a further decrease in melting temperature, and increase 
in critical temperature over the Lennard- Jones case. No new 
phase behaviour is observed. 



0.292 to 0.500 in steps of 0.042 locates the melting tem- 
perature at T* = 0.315 ± 0.01. Again the uncertainty 
is dominated by statistical error in the liquid reference 
free energy. Gibbs-Duhem integration initiated from this 
point proceeded in the direction of positive temperature, 
producing the phase diagram presented in Fig. [7| 



C. Phase Behaviour for A — 0.55e 

It is clear that the choice of A = e/2 has not captured 
the interesting phase behaviour expected in this region, 
specifically the transition from sh to fee structure. In- 
creasing A to 0.55e widens the energy difference between 
these two structures and should therefore manifest the 
transition at higher temperature and pressure. A study 
of the condensed phase behaviour of this model therefore 
seems appropriate. 

To locate the pressure of the sh-fee transition at zero 
temperature, the optimised enthalpy as a function of 
pressure was plotted for both structures. This enthalpy 
was obtained by CG minimisation with respect to both 
atomic positions and cell vectors. The two enthalpy 
curves intersect at approximately P* = 8.9. At this pres- 
sure the sh structure is of lower density. By simple con- 
sideration of the Clausius-Clapeyron equation these re- 
sults require that the transition at higher temperature oc- 
curs at lower pressure. This provides a range over which 
the finite temperature transition can be sought with free 
energy calculations. A finite temperature point on the 
sh-fee phase boundary was sought along the T* = 0.167 
isotherm. It should be stressed that density responses 
to pressure and temperature in the sh structure occur 
anisotropically. Expansion along each crystal direction 
must be considered separately when constructing a cell 
for thermodynamic integration. Free energies were cal- 
culated for pressure in the range P* — 0.9 to 1.3. Ex- 
trapolation between sampled points yields a transition at 
P* = 1.15 ±0.01. This error is estimated using the devi- 
ation from ideal finite-size scaling of the leading term in 
1/N when computing finite-size corrections. 

Melting of the sh structure has been studied along the 
P* = 0.118 isobar, conducting free energy calculations 
close to the A = e/2 melting line. A reference point for 
computing liquid free energies was taken at T* — 1.171 
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with chemical potential /i* = —9.474 and Helmholtz free 
energy /* = —12.95 ± 0.01 per atom. The intersection 
of the sh-solid and liquid free energies provides a melt- 
ing temperature of T* = 0.343 ± 0.008. By a similar 
process, the fee melting temperature at P* = 2.373 was 
determined as T* = 0.516 ± 0.008. 

These two melting points, plus the identified sh-fee 
transition pressure at T* — 0.167 provide starting points 
from which to begin Gibbs-Duhem integration. The three 
Gibbs-Duhem series are shown in Fig. |H1 The sh-fee 
meets the intersection of the two melting lines. The 
same triple point is located by the intersection of any 
two series, and is independently confirmed by the third 
to within a small error. The location of the triple point 
is hence T t * p = 0.358 ± 0.002, P* p = 0.54 ± 0.02. 

The Gibbs-Duhem information also confirms that both 
melting points measured are thermodynamically stable, 
information which was not available from the above free 
energy calculations alone. 
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FIG. 9: Optimised enthalpy per atom for fee, sh and sc struc- 
tures as a function of pressure at zero temperature in the 
A — e case. Both sc-fee (metastable) and sh-fee phase transi- 
tions are indicated. 



D. Phase Behaviour for A = e 

Here we expect the sh structure to dominate at low 
pressure. The transition to fee is expected to occur at 
significantly higher pressures than the A = 0.55e case. 

A suitable bimodal number density histogram from 
which to begin the histogram re- weighting/multi- 
canonical sampling procedure was identified at T* = 
1.671 when using a chemical potential of \i* =-14. 137. 
Density data produced from stepping along the coexis- 
tence curve identifies the critical temperature as T* — 
1.680 ± 0.004 when extrapolated to zero density differ- 
ence. This leads to a critical density of p* = 0.2716 ± 
0.0005. Below temperatures of T* = 1.450 the coexis- 
tence curve was traced with Gibbs-Duhem integration. 

Zero pressure free energy calculations for the relevant 
solid structures have been reported in Ref. and will 
not be repeated here. These confirm that the sh struc- 
ture is thermodynamically stable at low pressure. Zero 
temperature phase transitions have been located by CG 
enthalpy optimisation of the relevant structures at var- 
ious pressures. Plots of enthalpy against pressure for 
the sc, sh and fee structures are shown in Fig. The 
metastable sc-fee transition (which is metastable) occurs 
at a pressure of P* = 3.63, with the sh-fee transition at 
P* = 13.31. The sh-fee transition at finite temperature 
was located by performing six Einstein crystal calcula- 
tions for each phase in the pressure range P* = 10 to 17 
at a temperature of T* = 0.292. The metastable sc-fee 
transition has not been explored at finite temperature. 

As noted in section IIIII the NPH two-phase method 
is useful in locating the sh melting temperature. These 
simulations in the A = e case have been described in 
Ref. 23], leading to a melting temperature of T* — 
0.609 ± 0.002 at P* = 0.237. To confirm this with 
free energy calculation, a liquid reference point at T* — 
1.500 with /i* = —12.368 was taken. Grand canoni- 



cal Monte-Carlo simulation reveals (N) — 366.1 ± 0.5, 
and (P*) = 3.59 ± 0.01 under these conditions. The 
Helmholtz free energy of the reference point is hence com- 
puted as /* = -16.99 ± 0.01 per atom. 

Based on the information provided by the two-phase 
simulations, the sh melting temperature was sought along 
the P* = 0.237 isobar. Free energies for both phases were 
computed at points between T* = 0.471 and 0.671 in 
steps of 0.041. Interpolation to the intersection provides 
a melting temperature of T* = 0.614 ± 0.009 which is 
in agreement with the two-phase NPH result. The fee 
melting has been studied at a pressure of P* = 24.86 by 
a similar procedure, locating the melting temperature at 
T* = 1.851 ±0.009. 

As with the A — 0.55e potential, the three pairs of free 
energy calculations have been used as starting points for 
Gibbs-Duhem integration. Any two of the resulting three 
series locate the same triple point to within the error of 
the initial free energy calculations. The triple point is 
located at T* p = 1.04 ± 0.01, P t * p = 11.7 ± 0.2. The 
sh melting curve can be traced deep into the fee region 
with a large range of metastability. Within this range, a 
maximum in the melting curve appears at slightly higher 
pressures than the sh-fee triple point. 

The full phase diagram for this potential is plotted in 

Fig. nni 

E. Phase Behaviour for A = 3e/2 

In this potential the sh melting curve may exhibit a 
maximum in the stable regime. In addition, an isostruc- 
tural fcc-fcc transition is expected at positive pressure. 

Rather than employing the expensive process of tracing 
fluid isotherms, a starting point for the multi-canonical 
sampling procedure was obtained by extrapolation from 
previous A values. Fig. ^] shows a plot of critical 
temperature against the outer Gaussian well depth. A 
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FIG. 10: Phase diagram of the A = e potential (inset) 
in both the pressure-temperature and temperature-density 
planes. The sh phase now dominates at low pressure. A 
maximum in the sh melting curve is marginally preempted by 
the fee phase. A further increase in critical temperature with 
A is observed. 



quadratic fit to the previously identified critical temper- 
atures yields T* — 2.142 when extrapolated to A = 1.5e. 
Based on this estimate, a suitable bimodal histogram 
from which to begin tracing the liquid vapour coexistence 
curve was identified at T* = 2.083 with a chemical po- 
tential of fj,* = —17.705. Extrapolation of the resulting 
multi-canonical density data identifies the critical point 
at T* = 2.12 ± 0.04, p* = 0.265 ± 0.002. 

The isostructural transition was first located at zero 
temperature using CG enthalpy minimisation, optimising 
both phases at pressures between P* = 0.237 and P* = 
2.37 in steps of 0.12. The resulting enthalpy per atom 
for both phases is shown in Fig. ^| The intersection 
reveals a transition pressure of P* = 1.09. Note that 
above P* = 1.4 the lower density phase collapses to the 
higher density structure during optimisation, indicating 
mechanical instability. 

To locate the transition at finite temperature, free en- 
ergy calculations were performed for both phases along 
the T* — 0.208 isotherm. The chemical potentials de- 
rived from these free energies are plotted in Fig. H3 
locating the transition pressure for this temperature as 
P* = 1.27. As with sh-fee transitions the error on this 
value is small, being approximately 0.01. Both free en- 
ergy series are corrected for finite-size effects. 

The isostructural transition predicted by IScala et alJ 
[l3| has now been identified as fcc-fcc at low tempera- 
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FIG. 11: Critical behaviour for the A — 3e/2 potential. 
The critical temperatures identified for smaller values of A 
are shown above. A quadratic extrapolation to the current 
A is shown. Histograms resulting from the subsequent re- 
weighting and multi-canonical sampling procedure are shown 
below. 



ture. To locate the extent of the fcc-fcc transition, Gibbs- 
Duhem integration has been employed in the direction of 
increasing temperature. 

The resulting series indicates an increase in transition 
pressure with increasing temperature. However, after 
just six steps the lower density fee phase becomes me- 
chanically unstable. At this point a large density differ- 
ence between the two phases still exists. The isostruc- 
tural transition does not end in a critical point, but ter- 
minates at the low density fee spinodal line. 

The sh-fee transition was located by the same proce- 
dure as in section UVDl The zero temperature transition 
is located by CG enthalpy minimisation at P* — 25.27. 
The transition along the T* = 0.21 isotherm was then 
determined by computing the free energy of both phases 
at ten pressures between P* = 23.67 and P* = 26.04. 
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FIG. 12: Metastable isostructural transition in the A = 3e/2 
potential. The enthalpy per atom at zero temperature is plot- 
ted along with the results of free energy calculations along the 
T* = 0.21 isotherm. 

The transition is located at P* = 25.17. The estimated 
error in this pressure is of order 0.01. 

The sh and fee melting temperature were sought along 
the P* — 0.237 and 23.67 isobars respectively. No two- 
phase simulations were performed in this case. A suitable 
range over which to perform free energy calculations was 
estimated from trends in earlier data. A reference point 
at T* — 1.500 with fi* = —12.368 was used in computing 
liquid free energies along both isobars. Grand canoni- 
cal Monte-Carlo simulation provides (N) — 415.6 ± 0.5, 
and (P*) = 4.52 ± 0.01, leading to a Helmholtz free 
energy of /* = —20.79 ± 0.04 per atom. Free energy 
plots along each isobar locate the sh melting tempera- 
ture as T* = 0.80, with the fee melting at T* = 1.45. In 
both cases the uncertainty is approximately 0.01. Gibbs- 
Duhem series for the sh-fee, sh-liquid and fcc-liquid tran- 
sitions have been computed to trace the remainder of the 
phase diagram. As can be seen in Fig. the sh melt- 
ing curve maximum now lies within the stable regime. 
The fee melting curve intersects at slightly higher tem- 
perature. There is hence a small range visible in both 
the temperature-pressure and temperature-density pro- 
jections for which melting is reentrant. The metastable 
fcc-fcc transition is also plotted in Fig. ^| The termina- 
tion of this lies well below the melting temperature. 



V. LIQUID ANOMALIES 

In mapping the above phase diagrams, no evidence of 
a thermodynamically stable liquid- liquid phase transition 
has emerged. In particular data from Gibbs-Duhcm inte- 
gration reveal that density is continuous along the liquid 
side of all melting and vaporisation curves. In addition 
no third fluid peak has emerged during multi-canonical 
sampling. A stable LLPT must meet either the melting 
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FIG. 13: Phase diagram of the A = 3e/2 potential (in- 
set). Pressure-temperature and temperature-density projec- 
tions are shown. The maximum in the sh melting curve 
is clearly visible. The dashed line indicates the metastable 
isostructural transition. An increase in both melting and crit- 
ical temperatures is observed over the A = e case. 



or vaporisation curves at a triple point and can hence be 
ruled out. 

The possibility of liquid anomalies has been investi- 
gated with a considerable number of constant pressure 
Langevin dynamics simulations. Particular attention has 
been focused on the A — 3e/2 liquid, in the region where 
the solid is less dense than the liquid. A total of 74 simu- 
lations have been conducted over the range P* = 16.5 to 
P* = 26.0, T* = 1.42 to 1.71. Each simulation is of dura- 
tion t* = 68 at equilibrium with a system size N = 500. 
No anomalies are observed in the density (Fig. 114(1 . heat 
capacity, bulk modulus or diffusion coefficient. This is to 
be expected. In contrast to the two-dimensional case in 
ref [lg > the melting transition is highly first-order. Nu- 
cleation of solid clusters is therefore strongly inhibited 
by finite-size effects. Accelerated sampling methods may 
prove useful in this region. It should be noted that the 
isostructural transition in the A = 3e/2 potential termi- 
nates at the ld-fcc spinodal line below the glass tempera- 
ture. This has been located using the method in ref [37| 
as T* ~ 0.47 at P* — 1.5. No influence on the super- 
cooled liquid is therefore expected. This is confirmed by 
a sequence of 12, N = 500 simulations (also of duration 
t* = 68) along the T* = 0.44 and T* = 0.56 isotherms. 
No anomalies or evidence of a LLPT are found in the 
metastable liquid. 
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FIG. 14: Liquid density close to the sh melting line for A — 
3e/2. Density against temperature along several isobars is 
shown with accompanying error bars. 



VI. LIQUID STRUCTURE 

Finally, we briefly examine the structure of the liquid.. 
The pair correlation function g(r) was computed under 
a variety of conditions and studied as a function of the 
parameter A. Results at T* = 1.6 are shown in figure IT51 
for densities p* —0.7 and p* — 1.0. These were computed 
from canonical ensemble molecular dynamics simulations 
using 500 atoms over a duration of t* — 170. 

At the lower density, increasing the depth of the Gaus- 
sian well leads to the emergence of a peak not present in 
the Lennard- Jones (A — 0.0) case. This corresponds to 
an coordination shell lying at an intermediate radial dis- 
tance between the first and second Lennard- Jones shells. 
The occupation of this shell occurs at the expense of the 
innermost peak. The effect of increasing the A parame- 
ter is hence to shift a significant fraction of first-nearest 
neighbors into the intermediate coordination shell. The 
second Lennard-Jones coordination shell is seemingly un- 
affected. 

At the higher density, the first coordination shell is 
unchanged on increasing A. The intermediate peak is less 
pronounced and emerges at the expense of narrowing the 
second Lennard-Jones coordination shell. Simulations at 
intermediate densities indicate that although significant, 
the change in liquid structure occurs continuously and 
does not correspond to a phase transition. 



VII. CONCLUSIONS 

We have mapped the liquid-gas, melting and solid-solid 
transitions in a sequence of core-softened pair potentials 
with increasing outer well strengths. The liquid-gas criti- 
cal temperature has been seen to increase approximately 
quadratically with increasing A. The pressure of the crit- 
ical point also increases. 

As A is initially increased, the fee melting temperature 
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FIG. 15: Evolution of the pair correlation function g(r) with 
A at two densities. In each case the temperature is T* — 1.6. 



is found to reduce, perhaps contrary to the expected be- 
haviour. This decrease in melting temperature can be 
understood by examining Fig. ^ For A values in this 
region, the effect of the Gaussian is to widen the existing 
Lennard-Jones minimum. This allows larger fluctuations 
about equilibrium lattice positions for a given tempera- 
ture. The well known empirical rule of Lindemann |3S| 
states that melting will occur when the r.m.s. fluctuation 
is ~ 15% of the nearest-neighbour distance. This will oc- 
cur at lower temperatures for wider potential wells. 

As the Gaussian outer minimum becomes distinct from 
the Lennard-Jones minimum, the simple hexagonal struc- 
ture becomes lower in energy due to a large number of 
third-nearest neighbours close to rrj. A transition to fee 
under pressure has been observed and is seen to intersect 
the melting line. The resulting simple hexagonal melting 
temperature increases with increasing A. The pressure 
of the sh-fee transition increases with increasing A, as 
expected from the larger enthalpy difference. 

We have also seen that the sh structure exhibits a max- 
imum melting temperature, which for larger A values is 
manifested in the thermodynamically stable regime. The 
predicted isostructural transition has been observed for 
the fee structure only. Energy minima are in fact also 
observed for two volumes when imposing simple cubic 
symmetry. The higher energy structure is however me- 
chanically unstable when atoms are permitted to move 
away from lattice sites. The fcc-fcc transition does not 
approach the melting line. 

The low pressure phase diagram and sublimation have 
not been studied here. The phase behaviour is expected 
to be unchanged from the Lennard-Jones case in this 
region. The hexagonal close-packed structure has also 
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not been studied in detail. For the range of parame- 
ters studied the fcc-fcc transition is not resolvable with 
the methods used here. The solid should be considered 
'close-packed' in regions where the fee structure has been 
calculated as energetically favourable. We note that the 
lattice-switch Monte-Carlo method [sti Ecj has recently 
been employed to resolve the fec-hep transition in the 
Lennard-Jones potential |4£j an d could be used to simi- 
lar effect in these systems. 

The two dominant sources of error in this study are 
the finite-size error in the solid and the statistical error 
associated with computing a liquid reference point for 
thermodynamic integration. The former has been largely 
corrected for by repeat calculations with larger system 
sizes. The latter dominates over the liquid finite-size er- 
ror, but has been controlled to an acceptably low level. 
All phase boundaries shown can be considered accurate 
to within temperatures of AT* = ±0.05 and pressures of 
AP* = ±0.1. 

The failure of the isostructural fcc-fcc transition to ex- 
tend into the supercooled liquid suggests that potentials 



of this specific function form (Lennard-Jones plus Gaus- 
sian) are not capable of reproducing the unusual phase 
behaviour of water. Variation of the outer well posi- 
tion and width may yield significantly different phase be- 
haviour to that reported here. The range of parameters 
over which the energetic ordering of phases is unchanged 
is currently under investigation. The temperature differ- 
ence between the ld-fcc spinodal and the glass transition 
is small. It may therefore be possible to extend the fcc- 
fcc transition beyond the glass temperature using small 
changes within this parameter range. 
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